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Abstract 

In this paper we propose on continuous level a class of domain decom- 
position methods of Robin-Robin type to solve the problems of unilat- 
eral contact between elastic bodies with nonlinear Winkler covers. These 
methods are based on abstract nonstationary iterative algorithms for non- 
linear variational equations in reflexive Banach spaces. We also provide 
numerical investigations of obtained methods using finite element approx- 
imations. 
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1 Introduction 

Thin covers from another material are often applied in engineering to improve 
the functional properties of the surfaces of components of machines and struc- 
tures. On the other hand, thin covers with certain mechanical properties are 
used for modeling of real microstructure of the surfaces, adhesion and glue 
bondings [rl [13. Ha]. 

The classical methods for solution of contact problems for bodies with thin 
covers are grounded on integral equations and are reviewed in work . Nowa- 
days, one of the most effective numerical methods for such contact problems are 
methods, based on variational formulations and finite element approximations. 
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Efficient approach for solution of multibody contact problems is the use of 
domain decomposition methods (DDMs). Many DDMs for contact problems 
without covers are obtained on discrete level 

HE!. Amon g DDMs, proposed 
on continuous level for contact problems without covers are methods presented 
in P, 0, [l2| ■ Domain decomposition methods for solution of problem of ideal 
contact between two bodies, connected through nonlinear Winkler layer are 
proposed in d • These methods are based on saddle-point formulation and 
conjugate gradient methods. 

In current contribution we consider the problem of unilateral contact be- 
tween bodies with nonlinear Winkler covers. We give variational formulations 
of this problem in the form of nonlinear variational inequality on convex set and 
variational equation in the whole space, and present theorems about existence 
and uniqueness of their solution. Furthermore, we propose on continuous level a 
class of parallel domain decomposition methods for solving the nonlinear varia- 
tional equation, which corresponds to original contact problem. In each iteration 
of these methods we have to solve in a parallel way linear variational equations 
in separate bodies, which are equivalent in a weak sense to linear elasticity prob- 
lems with Robin boundary conditions on possible contact areas. These DDMs 
are based on abstract nonstationary iterative methods for variational equations 
in Banach spaces. They are the generalization of domain decomposition meth- 
ods, proposed by us earlier in 0, M [Toj j for unilateral contact problems without 
covers. Some particular cases of proposed DDMs can be viewed as a modifi- 
cation of semismooth Newton method [7\ . The numerical analysis of obtained 
DDMs is made for plane contact problems using finite element approximations. 



2 Statement of the problem 

Consider a unilateral contact of N elastic bodies £l a C K 3 with piecewise smooth 
boundaries r Q , a — l,2,...,N (Fig. la). Suppose that across each contact 
surface there is a nonlinear Winkler layer. Denote 57 = (J^Li 




2h 

p 2 ,=o 

p 22 =0 

h 

P„=0 
P, 2 =0 



b) 







' P 2 ,=0 p 22 =-q 




IHUUli 




n 2 


u 21 =0 




p 22 =0 


V b' 


u«=0 




P,2=0 


[tttttttt' 





Pi,=0 p 12 =q 



Figure 1: Unilateral contact between several elastic bodies through nonlinear 
Winkler layers 

A stress-strain state in point x = (x\, X2, x$) T of each solid il a is described 
by the displacement vector u a — u a j e.; , the tensor of strains e a = e a ij ej Bj 
and the tensor of stresses a a = <r a ij e.; e 3 . These quantities satisfy the following 
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relations: 

dx 



]l ;; j(x) +/ai(x)=o, x£^ Q , 1=1,2,3, a) 

j=i j 



o- Qij -(x) = ^ C ai jfej(x)e afeJ (x), xG^«, i,j = l,2,3, (2) 



fe,i=i 



1 / du a i (x) 9u Q ,• (x) \ ^ „ „ „ 



(3) 



where / a , are the components of volume forces vector f Q = f a i&i, and C a ijki 
are symmetric elasticity constants, which are bounded in the following sense: 
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(3b a ,c a > 0) (Vx) I b a £ lij ^ X! Ca ^ kl <c a J2 ) ■ ( 4 ) 

[ i,j=l i,j,k,l=l k,l=l J 

Introduce on boundary T a a local orthonormal coordinate system £ a , r) a , n a , 
where n Q is an outer unit normal. Then the vectors of displacements and stresses 
on r Q can be written in the following way: u a — u a j £ a +u an rj a +u 

an Hqi & ol — 

o a ■ n a = cr Q? £ Q + cr av r) a + a an n a . 

Suppose, that the boundary T a consists of three disjoint parts: 
r Q - KUKiJSa, K = T", T u a ± 0, S a ^ 0. On the part V" homogenous 
Dirichlet boundary conditions are prescribed, and on the part T^ we consider 
Neumann boundary conditions: 

u«(x) = o,xer: ; <x Q (x) = Pa (x), xer». (5) 

The part S a — {Jp eBa S a p, C\fjeB a = is the possible contact area of 
body Sl a with the other bodies. Here S a p is the possible unilateral contact area 
of body Sl a with body Sip, and B a C {1,2, N} is the set of the indices of 
all bodies in contact with body Sl a . We assume that the surfaces S a p C T a 
and Sp a C Tp are sufficiently close (S a p ~ Sp a ), and n Q (x) w — n ( g(x'), 
x G S a p, x' = P(x) G Sp a , where P(x) is the projection of point x on S a p. Let 

d a p(x) = ±||x — x'|| = ±\Jj2i=i — x 'i) 2 be a distance between bodies Sl a 
and Sip before the deformation. 

We suppose that possible contact areas S a p and Sp a , P G B a , a — 1,...,N 
have nonlinear Winkler covers. Total compression w a p of these covers is related 
with normal contact stress as follows: <7 an (x) = a 0n (x') = g a p {w a p(x}), x G 
S a p, x' G Sp a , where g a p is given nonlinear continuous function, which satisfy 
the next conditions: 

g a p(0) = 0, (Vy,z) {y < z => g a p(y) < g a p{z)} , (6) 

{3M a p >0) (Vy,z) {\g a p(y)- g a p(z)\ < M a p \y - z\ } . (7) 

On possible contact zones S a p, G B a , a = 1,2, ...,7V we consider the 
following unilateral contact conditions through nonlinear Winkler layers: 

(7 q4 (x) = CT^(x') = , <T an (x) = <Tp n (x.') = 0, (8) 

CTan(x) = C7 i a n (x') = a/3 (tU aj8 (x)) < , (9) 
Man(x) + U/3n(x') + «J a/ 3(x) < d Q(3 (x) , (10) 

[u a „(x) + w 1 g n (x') + w a p{x) - d a p(x)](T an (x) = 0, x' = P(x), x G S Q/3 . (11) 
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3 Variational formulations 

For each body £l a consider Sobolev space V a — [i/ 1 (f2 Q ,)] 3 and the closed sub- 
space V Q ° = {u a e V a : u Q = on r^}. All values of the elements from these 
spaces on the parts of boundary T a should be understood as traces. The trace 
of element u Q e V a on the part should belong to space [7? 1 / 2 (r^)] 3 , and 
the trace of element from V® on the part S Q = int (r Q \ T^) should belong to 

[< /2 (S a )] 3 . 

Define Hilbert space Vb = Ila=i ^ a w ith scalar product 
( u ' v )v = Ea=i ( u a,v a ) Va and norm ||u|| Vo = (u,u)^ 2 , u, v G V . Moreover, 
introduce following spaces W = ]J {a ^ }gQ H^ 2 (E a ) = {w = (w a/3 )J a ^ }eQ : 

w a p G Hqq 2 } and U = Vq x W = {U = (u,w) T : u G Vb, w G W}, where 
Q = {{a,(3}: a G {1, 2, N} , (3 e B a }. 

In space Uq consider the closed convex set of all displacements, which satisfy 
nonpenentration contact conditions: 

K = {U 6 Uq : u an + up n + w a p < d a p on S a p, {a, (3} G Q }, (12) 

where w Q „ = n Q • u a G i? ^ 2 (E; Q ), w a p,d al3 G H^ 2 (E a ). 

Let us introduce bilinear form A(u, v), such that A(u, u) represents the 
total elastic deformation energy of the bodies, linear form L(u), which is equal 
to external forces work, and nonquadratic functional £f(w), which represents 
the total deformation energy of nonlinear Winkler layers: 



A(u,v) = y~] a Q (u Q ,v Q ), a Q (u Q ,v a )= / a a (u a ) :e a (v a )dCl, (13) 
L (u) = V] Z Q (u Q ), l a (u a )= t a -u a dQ+ p a -u a dS , (14) 

a = l ^ r a 

#( w ) = ^2 / g a (s(z)dz 



{a,/3}£Q' 



dS, u,veVo, weW, (15) 



where f Q G [L 2 (n a )} 3 , p Q G [H 00 1/2 (S Q )] 3 , a = l,2,...,iV. 

We have shown, that if condition (HJ holds then bilinear form A is sym- 
metric, continuous and coercive, and nonquadratic functional H is Gateaux 
diffcrentiable: 

H'(w,z)= ^2 / 9af3{w a p) z a/3 dS , w,zef . (16) 

Theorem 1. Suppose that conditions Op, (0), (0) hold. Then problem (Q|)- 
(0), J3J), {SP-fliP ^ as an alternative weak formulation as the following mini- 
mization problem: 

F(U) = A (u, u)/2 - L (u) + H(w) -> nun . (17) 

Moreover, there exists a unique solution of problem |J7[ ), and i/izs problem is 
equivalent to the following nonlinear variational inequality on set K : 

F'(U, V - U) = A (u, v - u) - L (v - u) + #'(w, z - w) > , V (v, z) T G X . 

(18) 
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Except this variational formulation, we also have proposed another weak 
formulation of original contact problem in the form of nonlinear variational 
equation. 

Let us introduce the following nonquadratic functional in space Vq: 



{a,/3}eQ' 



J ( u )= Yl / / 9^p{ z )dz 

r_TT^^JS„« JO 



dS, u e V , (19) 



where g a p(z) = {0, z > } V { g a p{z) , z < } is nonlinear function. 
Functional J(u) is nonnegative and Gateaux differentiable in Vq: 

J'(u,v) = - ^2 / 9ap{d a p - u an - up n ) [v an +vp n ] dS . (20) 

{a,0}EQ Js °>? 

We have shown that if conditions © and ([7]) hold, then Gateaux differential 
J'(u, v) satisfies the following properties: 

(VuGVb) (3R>0) (VveVo) {|^(u,v)| <fl||v|| Vo }, (21) 

(3£>>0)(Vu,v ) weVb){|J / (u + w ) v)-J / (u,v)|< J D[[v|[ vb ||w|| Vb } ! (22) 

(Vu,v G Vb) { J'(u + v,v) - J'(u,v) > 0}. (23) 

These properties helped us to prove the next theorem. 

Theorem 2. Suppose that conditions (0J) and ^ hold. Then the 
contact problem (QP-lEJ), ([SJ), is equivalent to problem 0), (0) 

wf/j the following nonlinear boundary value conditions on the possible contact 
areas: 

<r an (x) = cr/3„(x') = g~p {d a p(x) - u an (x) - u Pn (x')) , x' = P(x) , x G S a p , 

(24) 

and it is equivalent in weak sense to the next nonquadratic minimization prob- 
lem: 

Fi (u) = A (u, u)/2 - L (u) + J(u) min . (25) 



Moreover, problem I125\) has a unique solution and is equivalent to the next 
nonlinear variational equation in space Vq: 

F 1 (u,v) = A(u,v) + J'(u,v)-L(v)=0, VveVb, ueVb. (26) 



4 Nonstationary iterative methods 

In reflexive Banach space V consider an abstract nonlinear variational equation 

$(u,v)=y(v), VveV, UGF, (27) 

where $ : V x V — ?> R is a functional, which is linear in v, but nonlinear in 
u, and Y : V — > M is linear continuous form. For numerical solution of (|2"T)l 
consider the next nonstationary iterative method 0, fllj j : 

G fc (u fc+1 ,v) = G fe (u fc ,v) - 7 fc [$(u fe ,v) -y(v)] , fc = 0,l,... , (28) 
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where G k : V x V — > K are some given bilinear forms, ]' £ I are iterative 
parameters, and u fe 6 V is the fc-th approximation to the exact solution of 
problem (|2"7|) . 

Theorem 3. [B[ Suppose that functional <£> satisfies the following properties: 

(Vu G V) (3i?$ > 0) (Vv G V){ |*(u, v)| < i?$ ||v|| v } , (29) 

(3I>*>Q)(Vu > v > weV){|$(u + w > v) -$(u,v)| < I>*||v|| v ||w|| v }, (30) 

(3fl$>0)(Vu,v6V) {*(u + v,v)-$(u,v)>B«||v||y}. (31) 

Then nonlinear variational equation \2ty has a unique solution u £ V. In 
addition, suppose that bilinear forms G k , fc = 0, 1, ... are symmetric, continuous 
with constant Mq > 0, coercive with constant Bq > 0, and the next conditions 
hold: 

(3fc £ N ) (Vfc > fc ) (Vu e V) { G fe (u, u) > G fe+1 (u, u)} , (32) 
(Be G (0, 7 *), 7 * = B^B* G /D%) (Bfc) (Vfc > fc x ) { 7 fe G [e, 2 7 * - e] } , (33) 
where {u fe } C is obtained by iterative method 128]) . 

5 Domain decomposition schemes 

Now let us apply nonstationary iterative method (1281) for solving nonlinear vari- 
ational equation (|26[) . which corresponds to original contact problem. This 
equation can be written in form ((27)) . where $(u,v) = A(u,v) + J'(u, v), 
y(v) = L (v) , u, v G V, V — Vo , and iterative method (|2"5|) applied to solve 
(J2HJ) rewrites as follows: 

G fc (u fc+1 ,v) = G fc (u fc ,v)- 7 fc [A(u fc ,v) + J'(u fc ,v) - L(v)] , fc - 0, 1,.... (34) 

Note, that in general case iterative method AMI) does not lead to domain 
decomposition. Let us propose such variants of this method, which involve the 
domain decomposition. At first, let us take bilinear forms G k in method ((34]) 
as follows: 

G fe (u,v) =d 2 Fi(u fc ,u,v) = A(u,v) +c» 2 J(u fc ,u,v), u,vGV , (35) 
<9 2 J(u fc ,u,v) = / X a /3 g'api.dafi-u^-u^) [u an + up n ] [v an + vp n ] dS, 

X a p = -[sgn(rf Q /j-uL-^™)]~ = {0, d a p-u k ar -u k n > 0}V{1, else}. (36) 

Here d 2 Fi (u k , u, v) and d 2 J(u k ,u,v) are the second subdifferentials of func- 
tionals F\ and J in point u fc G Vq. In the case when -f k = 1, fc = 0,1,..., 
iterative method (|3"4"|) with bilinear forms ([33)) corresponds to semismooth New- 
ton method for variational equation (f2l)|). However, this method does not lead 
to domain decomposition. 

Now, let us take bilinear forms G k in the following way: 

G k (u, v) = A (u, v) + X k (u, v) , u, v G Vq , (37) 
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^ fc (u,v) = ^ / ^f}9' a f3(d a p-u k an -u k 0n )u an v an dS, u,veVo, (38) 

where VS^( X ) ~ 1 x e I v { 0> x e ^ap\S k 8 } are characteristic functions 
of some given subsets S^g C S a p of possible contact areas. 

Let us show, that such choice of bilinear forms G k involves the domain 
decomposition. Introduce a notation u fc+1 = (u fc+1 — u k )/'y k + u k £ Vq. Then 
iterative method (I34[) with bilinear forms (|37l) can be written in such way: 

A(u k+1 ,v)+X k (u k+1 ,v) = L{v)+X k (u k ,v)- J'(u fc ,v), VveVb. (39) 

u fc+1 = 7 fc u fc+1 + (1 - 7 fe ) u fc , k = 0,l,.... (40) 

Since the common quantities of the subdomains are known from the pre- 
vious iteration, variational equation (1391) splits into N separate equations in 
subdomains f2 Q , and iterative method (131)1) - (|4T)|) can be written in the following 
equivalent form: 

a a {u k+ \v a )+ / 4>Zf) g'cpidaf} ~u k n -u k p n )u 1 £ l 1 v an dS = 

= l a (v a )+ ^ / i> k a g' aB (d aa -u k n -u k an )u k an v an dS + 

+ S / 9^{da0-u k an -u k Pn )v an dS, Vv a eV Q °, (41) 

< +1 =7 fe fi" +1 + (l-7 fc R, a = 1,2, ...,7V, fc = 0,l,.... (42) 

In each iteration k of method (14"T1) - (|4"2")) . we have to solve N linear varia- 
tional equations (1411) in parallel, which correspond to linear elasticity problems 
in separate bodies tt a with Robin boundary conditions on possible contact areas. 
Therefore, this method refers to parallel Robin-Robin type domain decomposi- 
tion schemes. 

By taking different characteristic functions ip k g, we can obtain different 
particular cases of domain decomposition method (|4"Tj) - (j4"2"|) . Thus, taking 
■0*3 (x) = (S k g = 0), Va,/3, Vfc, we get parallel Neumann-Neumann domain 
decomposition scheme. Other borderline case is when ip k B (x) = 1 {S k g = S a p), 
Va,/3, Vfc. 

Moreover, we can choose characteristic functions ip k 8 by formula (f31)|). i.e. 
ip k g = XaB- Numerical experiments, provided by us, have shown, that such 
DDM has higher convergence rate than other particular domain decomposition 
schemes. 

6 Numerical investigations 

Numerical investigations of proposed DDMs have been made for plane problem 
of unilateral contact between two isotropic bodies f2i and VL2, one of which has 
a groove (Fig. lb). The bodies are uniformly loaded by normal stress with 
intencity q — 10 MPa. Each body has length I = 4 cm and height h = 1 cm. 
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The Young's moduli and Poisson's ratios of the bodies are the same: E\ = 
E2 = 2.1 • 10 5 MPa, vi = V2 = 0.3. The distance between bodies is <ii2(x) = 



r{[l-(x 1 -l) 2 /b 2 } + } 3/2 , x G 5i2, where b = 1 cm, r = 5 
max{0, z}, S12 = {x = (xi 1 x2) T ■ x\ G [0, I], X2 = h}. 

Across possible contact area S12 there is a nonlinear Winkler layer. The 
relationship between normal contact stresses and displacements of this layer are 
described by the following power function: 

512 (wi2 (x)) = sgn (iui2(x)) Itoiafx)! 1 ' , x G S12 , where parameters B 

and a are taken from the intervals B G [ 10~ 6 cm/(MPa)°, 2 • 10~ 4 cm/(MPa) a ] , 
a G [0.1, 1]. For such choice of these parameters the nonlinear Winkler layer 
models a roughness of the possible contact surface [g| . 

This problem has been solved by DDM (|4"Tj) - (|4"2")) with stationary iterative 
parameters 7* = 7, Vfc and characteristic functions ipi 2 , taken by formula (1551) . 
i.e. ipi 2 = X127 F° r solving linear variational problems (|4Tj) in each iteration 
k we have used finite element method with 8192 linear triangular elements for 
each body. 

We have used the following initial guesses for displacements u? n (x) = 10~ 4 , 
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Figure 2: Relative error (a), and normal contact stress (b) 

At Fig. 2a the relative error p\ of displacement U2 n on different iterations 
k, obtained for B — 2.5 • 10 -5 cm/(MPa) a , a = 0.5, is represented for different 
values of parameter 7. Curves 1-9 correspond to 7 = 0.01, 0.02, 0.05, 0.6, 0.8 
(0.3), 0.9, 0.95, 0.98, 0.99. For these values of parameter 7, DDM (|?I ]) -(|i2 l) 
reaches the accuracy e u = 10~ 3 in 193, 124, 65, 7, 13, 27, 55, 134 iterations 
respectively. 

Thus, we conclude, that the best convergence rate reaches if 7 = 0.6. The 
convergence rate is good if 7 G [0.1, 0.9]. However, it becomes slow when 7 
is close to or to 1. For 7 = 0.98 the method is still convergent, but the 
convergence becomes nonmonotone. For 7 > 0.99 the method is not anymore 
convergent. 

At Fig. 2b the normal contact stress <j\ n = cr2„, obtained by DDM (|4lT) ((42)) 
for B = 10~ 5 cm/(MPa) a and different values of parameters a is represented. 
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Curves 1-4 correspond to numerical solution for a = 0.3, 0.6, 0.8, 1. Dashed 
curve represents the analytical solution, obtained in 13] for contact of two 
halfspaces without nonlinear layer. Here we conclude, that for small values of a 
(a < 0.3) the influence of nonlinear layer on the contact behavior is not so large 
and the numerical solutions are close to the solution without layer. However, 
for larger values of a (a > 0.5) the influence of nonlinear layer becomes more 
significant and can not be neglected. 

The positive feature of proposed domain decomposition methods are the 
simplicity of their algorithms. These methods have only one iteration loop, 
which deals with domain decomposition, nonlinearity of Winkler layers and 
contact conditions. 
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